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Abstract 



Charged particle production has been measured in Deep Inelastic Scattering (DIS) 
events using the ZEUS detector over a large range of Q 2 from 10 to 1280 GeV 2 . The 
evolution with Q of the charged multiplicity and scaled momentum has been investigated 
in the current fragmentation region of the Breit frame. The data are used to study QCD 
coherence effects in DIS and are compared with corresponding e + e~ data in order to test 
the universality of quark fragmentation. 
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1 Introduction 



A widely studied topic of perturbative Quantum Chromodynamics (pQCD) is the nature of the 
quark fragmentation process from the quark-antiquark pair created in e + e~ annihilation experi- 
ments. The purpose of this paper is to compare the fragmentation properties of the struck quark 
in Deep Inelastic Scattering (DIS) to those of the quarks produced in such e + e~ annihilation 
experiments, in order to test the universality of the quark fragmentation process. 

Operation of the HERA electron-proton collider has vastly extended the kinematic range for 
studies of the hadronic final state in DIS. The event kinematics of DIS are determined by the 
negative square of the four-momentum transfer, Q 2 = —q 2 , and the Bjorken scaling variable, 
x = Q 2 /2P-q, where P is the four-momentum of the proton. In the Quark Parton Model 
(QPM), the interacting quark from the proton carries the four-momentum xP. The variable 
y, the fractional energy transfer to the proton in its rest frame, is related to x and Q 2 by 
y = Q 2 /xs, where \fs is the electron-proton centre of mass energy. 

A natural frame to study the dynamics of the hadronic final state in DIS is the Breit frame 0. 
In this frame the exchanged virtual boson is completely space-like and is given by 
q = (0,0,0,— Q = —2xP B ) = (E,p x ,p y ,p z ) where P B is the proton momentum in the Breit 
frame. In the QPM the z-component of the momentum of the incoming quark is Q/2 before 
and —Q/2 after the interaction with the exchanged virtual boson. The particles produced in 
the DIS interaction can then be assigned to one of two regions. The current region corresponds 
to the direction of the outgoing struck quark in the QPM. In this paper we define the current 
region corresponding to particles with p z < 0, whilst the target region is defined by p z > 0. 

In e + e~ annihilation the two quarks are produced with equal and opposite momenta, ±-y/i/2. 
In the Breit frame of DIS, a quark is struck from within the proton with outgoing momentum 
—Q/2. In the direction of the struck quark the particle spectra are expected to have no 
dependence on x and a dependence on Q similar to those observed in e + e~ annihilation [Q, |3|. |J] 
at energy y/s = Q. 

The results from e + e~ annihilation support the inclusion of coherence effects in pQCD |5], || [7), |8[]. 
The phenomenon of coherence is a natural consequence of the quantum mechanical treatment of 
the parton cascade. As long wavelength gluons are unable to resolve individual colour charges 
of partons within the parton cascade, the available phase space for soft gluon emissions is 
reduced to an angular-ordered region, due to destructive interference. This leads to a number 
of important differences in the properties of the partonic final state relative to the incoherent 
case. The most notable of these are the slower rise in the multiplicity of partons with increasing 
energy and the modification of the logarithmic parton momentum spectra to an approximately 
Gaussian form, which is often referred to as the "hump-backed" plateau ||. 

Coherence effects are explicitly included in the Modified Leading Log Approximation (MLLA) 0] 
of pQCD. The MLLA calculations predict the parton multiplicity and the form of the momen- 
tum spectra from quark fragmentation at a given, sufficiently large, energy. The hypothesis of 
Local Parton Hadron Duality (LPHD) ||, which relates the observed hadron distributions to 
the calculated parton distributions via a constant of proportionality n ch , is used in conjunction 
with the predictions of the MLLA. Using this assumption the modifications due to coherence 
are therefore not only expected to be experimentally observed, but the MLLA calculations 
should also be directly applicable to data. 



1 



So far these effects have been studied at e + e~ annihilation experiments [10, 11, 12 1 using the 
average charged particle multiplicity, < n ch >, and the distribution of the scaled momentum, 
ln(l/x p ), where x p = 2p/y/s. The results from these experiments are in accordance with the 
MLLA and the assumptions of LPHD. 

Previous fixed-target DIS (up) measurements in the Breit frame have been performed in the 
range 1 < Q 2 < 45 GeV 2 with a mean v beam energy of 33 GeV [JL3|J . In the present paper 
the scaled momentum and charged multiplicity distributions of the hadronic final state are 
measured in the current region of the Breit frame as a function of x and Q in the range 
6 x 10~ 4 < x < 5 x 10~ 2 and 10 < Q 2 < 1280 GeV 2 with a mean energy of the electron beam 
measured in the proton rest frame ~ 10 4 GeV. Comparisons are made with Monte Carlo models, 
MLLA analytic calculations and e + e~ data. The DIS data were obtained in 1993 with the ZEUS 
detector at the HERA collider where electrons of energy E e = 26.7 GeV collided with protons 
of energy E p = 820 GeV, and correspond to an integrated luminosity of 0.55 pb _1 . 



2 The ZEUS detector and trigger 



ZEUS is a multipurpose magnetic detector which has been described elsewhere fl~4] , [15] . Here 
we give a brief description concentrating on those parts of the detector relevant for the present 
analysis. 

Charged particles are tracked by the inner tracking detectors which operate in a magnetic field 
of 1.43 T provided by a thin superconducting coil. Immediately surrounding the beampipe is 
the vertex detector (VXD) which consists of 120 radial cells, each with 12 sense wires (16| . 
The achieved resolution is 50 /im in the central region of a cell and 150 /mi near the edges. 
Surrounding the VXD is the central tracking detector (CTD) which consists of 72 cylindrical 
drift chamber layers, organised into 9 'superlayers' [0. These superlayers alternate between 
those with wires parallel (axial) to the collision axis and those inclined at a small angle to give 
a stereo view. The magnetic field is significantly inhomogeneous towards the ends of the CTD 
thus complicating the electron drift. With the present understanding of the chamber, a spatial 
resolution of 260 /zm has been achieved. The hit efficiency of the chamber is greater than 95%. 

In events with charged particle tracks, using the combined data from both chambers, recon- 
structed primary vertex position resolutions of 0.6 cm in the Z direction and 0.1 cm in the 
XY plane are measured. (The ZEUS coordinate system is defined as right handed with the 
Z axis pointing in the proton beam direction and the X axis horizontal pointing towards the 
centre of HERA. The polar angle 6 is defined with respect to the Z-direction.) The resolution 
in transverse momentum for full length tracks is <j Pt /pt = 0.005pr © 0.016 (for pt in GeV). 

The solenoid is surrounded by a high resolution uranium-scintillator calorimeter divided into 
three parts, forward (FCAL), barrel (BCAL) and rear (RCAL). Holes of 20 x 20 cm 2 in the 
centre of FCAL and RCAL are required to accommodate the HERA beam pipe. Each of the 
calorimeter parts is subdivided into towers which in turn are segmented longitudinally into 
electromagnetic (EMC) and hadronic (HAC) sections. These sections are further subdivided 
into cells, which are read out by two photomultiplier tubes. A detailed description of the 
calorimeter is given in 
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For measuring the luminosity as well as for tagging very small Q 2 processes, two lead-scintillator 
calorimeters are used . Bremsstrahlung photons emerging from the electron-proton interac- 
tion point (IP) at angles 7 < 0.5 mrad with respect to the electron beam axis hit the photon 
calorimeter at 107 m from the IP. Electrons emitted from the IP at scattering angles less than 
6 mrad and with energies, E' e , between 0.2E e and 0.9E e are deflected by beam magnets and hit 
the electron calorimeter placed 35 m from the IP. 

For events with the scattered electron detected in the calorimeter, the trigger was essentially 
independent of the DIS hadronic final state. The trigger acceptance was greater than 97% for 
Q 2 > 10 GeV 2 and independent of Q 2 pi. 



3 Event selection 



The offline selection of DIS events was similar to that described in our earlier publications 
[21-26]. Scattered electron candidates were selected by using the pattern of energy deposition 
in the calorimeter. The electron identification algorithm was tuned for purity rather than 
efficiency. In studies with Monte Carlo DIS events and test beam data the purity was estimated 
to be > 96 % for E' e > 10 GeV. 

The ZEUS detector is almost hermetic, allowing the kinematic variables x and Q 2 to be recon- 
structed in a variety of ways using combinations of electron and hadronic system energies and 
angles. Measurements made purely from the measurement of the scattered electron energy and 
angle, are denoted with a subscript e whilst those made purely from the hadronic system, by 
the Jacquet-Blondel method, are denoted by the subscript JB ^7|. The double angle method 



which calculates the kinematic variables from the scattered electron angle and the angle of the 
struck quark, 7^, as calculated from the final hadronic system are denoted by the subscript 
DA ||. 

In the kinematic region selected, the optimal method of determining Q 2 , x and the boost, (3, 
from the laboratory to the Breit frame is based on the DA method. 

For our final event selection we demanded 



• E' e > 10 GeV, to achieve a high purity sample of DIS events; 

• Q 2 DA > 10 GeV 2 ; 

• He < 0.95, to reduce the photoproduction background; 

• Vjb > 0.04, to give sufficient accuracy for DA reconstruction; 

• 5 = — cos$i) > 35 GeV, where the sum runs over all calorimeter cells. The 
measured energy of the cell is denoted by Ei and its polar angle with respect to the 
incident proton beam by 0j. Nominally, 5 should peak at twice the electron beam energy, 
53.4 GeV. This cut is used to remove photoproduction events and to control radiative 
corrections. 
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Furthermore we required 



• a primary vertex position, determined from VXD and CTD tracks, in the range 
—50 < Z vtx < 40 cm and a cut on the radial distance from the beamline 

Rvtx = \J X%t x + Yvtx < 10 cm; 

• the impact point (X, Y) of the scattered electron in the RCAL to lie outside a square of 
32 x 32 cm 2 centred on the beam axis, to ensure the electron is fully contained within the 
detector and its position can be reconstructed with sufficient accuracy; 

• no more than 5 GeV of energy deposition in the electron calorimeter of the luminosity 
detector, to exclude potential photoproduction background events. 



Following these cuts, the remaining photoproduction background was estimated to be ~1% 
using events tagged in the luminosity detector. The contamination from beam-gas background 
was estimated to be below 0.5% from unpaired electron and proton bunches. Finally, we rejected 
QED Compton scattering events and residual cosmic and beam-related muons. 

A total of 31.5K events was selected in this way corresponding to an integrated luminosity of 
0.55 pb _1 . Of these events about 8-10% p3f contain a large rapidity gap in the hadronic final 



state. These are characterised by rj max < 1.5, where r] max is the maximum pseudorapidity of 
any calorimeter cluster in the event. Here, the pseudorapidity is defined by rj = — ln(tan(0/2)) 
and a cluster is an isolated set of adjacent cells with summed energy above 400 MeV. 



4 Track reconstruction and selection 

Two programs for track finding and fitting have been developed independently which follow 
different strategies in the pattern recognition and track fitting stages. 

In the first approach, adopted for the final data analysis, the track finding algorithm starts 
with hits in the outermost axial superlayers of the CTD. As the trajectory is followed inwards 
to the beam axis, more hits from the axial wires and the VXD are incorporated. The resulting 
circle in the XY projection is used for the pattern recognition in the stereo superlayer pattern. 
The momentum vector is determined in a 5-parameter helix fit. 



The other track finding program is based on the Kalman filtering technique ||29|| . Seed tracks 
found in the outer layers of the CTD are extended inwards and points are added as wire layers 
of the CTD are crossed. The track parameters at each step are updated using the Kalman 
method. In the second step a Kalman fit to the points found in the pattern recognition phase 
is performed taking into account non-linear corrections to the measured drift time. Following 
the reconstructed CTD track inwards, CTD and VXD hits are associated to the track. The 
CTD tracks are merged to VXD track segments using the Kalman filtering algorithm. 

Multiple Coulomb scattering in the beampipe, and the walls of the VXD and CTD were taken 
into account in the evaluation of the covariance matrix. The vertex fit is performed with the 



fitted tracks using the perigee parameterisation p0[ . The vertex position is evaluated and the 
track parameters at the vertex are re-evaluated. 
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The reconstructed tracks used in this analysis are associated with the primary event vertex and 
have pr > 200 MeV and \rj\ < 1.5, corresponding to the polar angle region between 25° and 
155°. This is a region of good CTD acceptance where the detector response and systematics 
are best understood. For tracks within the fiducial volume defined by these cuts, the track 
reconstruction efficiency is ~ 95%. 

The Breit frame boost was reconstructed using the scattered electron, where the energy was 
determined using the DA formula and the polar and azimuthal angles were measured from the 
impact point on the calorimeter. The four-momentum vectors of the charged particles were 
boosted to the Breit frame, assuming the pion mass to determine the particle's energy, and 
were assigned to the current region if p z < in the Breit frame. 

Monte Carlo simulations were used to determine the acceptance for tracks in the current region 
as a function of (x, Q 2 ). The chosen analysis intervals in (x, Q 2 ) correspond to regions of high 
acceptance ^70% in the current region of the Breit frame. These intervals are commensurate 
with the resolution of x and Q 2 . 

The measured variables in the analysis are the inclusive mean charged multiplicity, <n c h>, in 
the current region of the Breit frame and the distributions of ln(l/x p ). Here, x v = 2p/Q, which 
is the momentum p of a track from the primary vertex measured in the Breit frame, scaled by 
Q/2, the maximum possible momentum (ignoring effects due to intrinsic kr of the quark within 
the proton). The peak position of the ln(l/x p ) distribution is denoted by hi{l/ x p ) max . 

Uncertainty in the reconstruction of the boost vector, j3, was found to be the most significant 
factor on the resolution of ln(l/x p ). Near the peak position of the ln(l/x p ) distributions the 
resolution is ~ 0.15 units, with no significant shift around the peak position, leading to a choice 
of bin width of 0.25 units. Migration of tracks from the current region to the target region was 
typically ~ 8%. The migration into the current region from the target fragmentation region is 
of a similar magnitude as the migration out. At low values of y, however, the level of migration 
is up to ~ 25%. In the low y region, the hadronic activity is low and the measurement of 7^ 
becomes distorted by noise in the calorimeter leading to a reduced x resolution and hence an 
uncertainty in f3. 



5 QCD models and event simulation 

Monte Carlo event simulation is used to correct for acceptance and resolution effects. The detec- 
tor simulation is based on the GEANT 3.13 |31]] program and incorporates our best knowledge 
of the apparatus. Details of the detector (MOZART) and trigger (ZGANA) simulation codes 
are given in 

Neutral current DIS events with Q 2 > 4 GeV 2 were generated using the HERACLES 4.4 
program [B2J which incorporates first order electroweak corrections. The Monte Carlo generator 
LEPTO 6.1 H, interfaced to HERACLES via the program DJANGO 6.0 |4j, was used to 



simulate QCD cascades and fragmentation. The parton cascade was modelled in different ways: 

• with the colour-dipole model including the boson-gluon fusion process (CDMBGF), using 
the ARIADNE 4.03 program. In this model coherence effects are implicitly included 
in the formalism of the parton cascade; and, 
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using the matrix element plus parton showers option (MEPS) within LEPTO, where 
coherence effects in the final state cascade are included by angular ordering of successive 
parton emissions. 



These models use the Lund string fragmentation model |36j for the hadronisation phase as 
implemented in JETSET 7.3 [37|. An additional sample of events was generated with the 
HERWIG 5.7 Monte Carlo 



where no electroweak radiative corrections were applied. The 
parton cascade includes coherence effects in a manner similar to that of LEPTO and a clustering 
model is used for the hadronisation [[7], 

For ARIADNE and HERWIG the parameterisation of the parton distribution functions was 
the MRSD'_ set ||40|| . The GRV parameterisation was used for the MEPS data set. These 
parameterisations have been shown to describe reasonably the HERA measurements of the 



proton structure function, F 2 [21, 42 1. 



Whilst these programs give a reasonable description of the observed energy flow J22[ they do 
not describe the excess of events observed with a large rapidity gap. The properties of these 
diffractive events are consistent with the exchange of a pomeron between the proton and virtual 



photon. The POMPYT Monte Carlo [43] models high-energy diffractive processes where the 
proton emits such a pomeron whose constituents take part in a hard scattering process with 
the virtual photon. A hard quark density distribution for the pomeron provides an acceptable 
description of this class of events p6| . The Monte Carlo program based on the model of 



Nikolaev-Zakharov [44] also gives an acceptable description of this class of events. 



For each of the above Monte Carlo programs the default parameters were used. The Monte 
Carlo event samples were passed through reconstruction and selection procedures identical to 
those for the data. 



6 Data correction 

The correction procedure is based on the detailed Monte Carlo simulation of the ZEUS detector 
with the event generators described in the previous section. Since the ARIADNE model gives 
the best overall description of our observed energy flow [^2| it is used for the standard corrections 
and unfolding of the distributions. 

The data are corrected for trigger and event selection cuts; event migration between (x, Q 2 ) 
intervals; QED radiative effects; track reconstruction efficiency; track selection cuts in px and 
rj; track migration between the current and target regions; and for the decay products of Kg 
and A decays which are assigned to the primary vertex. 

Correction factors were obtained from the Monte Carlo simulation by comparing the "true" 
generated distributions before the detector and trigger simulations with the "observed" distri- 
butions after these simulations followed by the same reconstruction, selection and analysis as 
the real data. The "true" distributions did not include the charged particle decay products 
of Kg and A and charged particles produced from weakly decaying particles with a lifetime 
> 10~ 8 s. Monte Carlo studies showed that up to 2% of the current region tracks from the 
reconstructed primary vertex are due to charged particles from the decay of Kg and A. In a 
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separate analysis, the shape of the Kg and A distributions are shown to be well reproduced by 
the Monte Carlo f45|. 



For the ln(l/x p ) distributions, correction factors were calculated for each ln(l/x p ) bin 

1 / dn \ I 1 / dn \ 



F(\n(l/x F 



N gen \d\n(l/x p )J gcn / A^ots \d\n(l/x PJ/ ohs 



where N gen (A^) is the number of generated (observed) Monte Carlo events in the (x, Q 2 ) 
interval. The overall correction factors are greater than unity and typically < 1.3 and are 
independent of ln(l/x p ) around the peak position. 

The correction procedure for the multiplicity distributions was performed in two stages. The 
dependence of this procedure on the Monte Carlo input is discussed in section [F[ The first stage 
was to correct for track reconstruction efficiency and track migrations between the current and 
target regions in each (x rec , Q 2 ec ) interval where x rec and Q 2 ec are the reconstructed values of x 
and Q 2 . In each interval a comparison was made between the observed multiplicity distribution 
P no { x rec, Qlec) i an d the generated distribution P np (x rec , Q 2 ec ). This comparison yielded the 
correction matrix M nptno (x rec ,Q 2 ec ) with elements defined by 

Q 2 ) °^ even ^ s w ^ n n p tracks generated when n Q tracks were observed 

n p ,n rec, r ec Q £ even |- s w [fa Uo tracks observed 

This matrix relates the observed to the generated distributions in each (x rec , Q 2 ec ) interval by 

Pn p {. x rec) Qrec) ~ -^n p ,n { X rec, Qrec) ' Pn (%rec, Qrec) ' 



The second stage corrected for migrations between the analysis intervals and for the acceptance 
of the event selection cuts using the correction factors 

q _ Ptruei^truei Qtrue) 
Prec{%rea Qrec) 

where the p's are normalised multiplicity distributions in each interval at the reconstructed 
level (rec) and at the generator level where no selection cuts are applied {true). 

The corrected multiplicity distribution was then calculated according to the formula 

Pn p (%truey Qtrue) ~ ^ ' M npjTlo (x rec , Q rec ) ' PnS^reci Qrec) 



The overall correction factors are in the range 1.25 — 1.5 and are dominated by the corrections 
from applying the matrix M n no (x rec , Q 2 ec ). In the lowest two Q 2 intervals the corrections from 
C are of the same size as those from the matrix. The corrected multiplicity distributions in 
each (x,Q 2 ) interval are henceforth denoted by P(n c h). 
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7 Systematic errors 



A number of systematic checks were performed in order to investigate the sensitivity of the 
corrected results to features of the analysis as described previously. These can be categorised 
under event selection, track and vertex reconstruction, track efficiency and selection, and data 
correction methods. 

Systematic uncertainties arising from remaining photoproduction background were studied by 
tightening the y e cut from 0.95 to 0.8 and the 5 cut from 35 GeV to 40 GeV. The effects of 
calorimeter noise on the reconstruction of x were checked by tightening the i/jb cut from 0.04 to 
0.05. The effect of a small mismatch in electron energy scale between the Monte Carlo and data 
was checked by scaling the energy of the calorimeter cells associated with the electron in the 
Monte Carlo. The overall systematic errors associated with the event selection were negligible. 

The analysis was repeated using the second approach to pattern recognition, track and vertex 
fitting (see section |) applied to both the Monte Carlo sample and the data. This was found 
to be a significant source of systematic error. The mean multiplicity results changed by 3 — 5% 
in general but the effect was up to 8%. Changes to ]n(l/x p ) max were typically 3% although the 
maximum shift observed was 9%. 

The efficiency of the tracking system was varied within the Monte Carlo by removing all hits 
from the vertex detector and 80% of hits from a random superlayer in the CTD. This is a 
conservative estimate of the systematic uncertainty due to inefficiencies within the detector. 
The largest change to the mean multiplicity was 3 — 5%. Removing the pt and r\ track selection 
cuts produced changes in the results which were within the statistical errors. (The intervals in 
(x, Q 2 ) were chosen in order to maximise the acceptance in the current region and therefore 
reduce the sensitivity to these cuts.) The changes due to varying the parameters in the vertex 
fit were again within statistical errors. 

The Monte Carlo-based correction techniques described in section |] are model dependent. In 
order to investigate the size of this effect the data were corrected with two independent Monte 
Carlo samples, MEPS and HERWIG. The differences in the simulation of the hadronic final 
state implicitly test the uncertainty due to the boost reconstruction which is derived from 
the DA method. The systematic uncertainties on Yr(1/ Xp) max and the mean multiplicity were 
typically between 4 — 7%. This is the largest single contribution to the overall systematic error. 
In the highest (x, Q 2 ) interval, where the statistics from Monte Carlo and data are limited, the 
largest change to the mean multiplicity was 17%. 

Inclusion of events with a large rapidity gap in the correction procedure was investigated by 
correcting the data with a sum of 90% ARIADNE and 10% POMPYT events which reasonably 
describes the DIS data p5] . Such events are a small contribution to the inclusive distribution 
and as a result the \rt(l/ x p ) max values are reduced by 1% whereas the mean multiplicities are 
more sensitive and increase by 3 — 5% due to the effects discussed in section [84]. 

The systematic deviations from the above checks were combined in quadrature to yield the 
quoted systematic errors. 



8 



8 Results 



8.1 Multiplicity 

The variables and distributions presented in the following sections have been corrected for 
detector and acceptance effects. The multiplicity distributions are binned in x and Q 2 as 
indicated in table [1], in the intervals 6 x 10~ 4 < x < 5 x 1CT 2 and 10 < Q 2 < 1280 GeV 2 , where 
the mean values of Q range from ~ 4 to 30 GeV. Figure || shows the measured multiplicity 
distributions in the current region of the Breit frame. With increasing Q the mean of the 
multiplicity increases. At low Q, there is a significant probability for no tracks to be found in 
the current region which we address later. The multiplicity distributions are compared to the 
ARIADNE Monte Carlo which reproduces the data in each (x, Q 2 ) interval. 

In order to compare the multiplicity distributions at various values of Q, the scaled multiplicity 
distributions, *&(z) =< n c h > -Pfoch), are plotted as a function of z = < ^ c ^ > in figure || |46[. 
For Q<7 GeV, the distributions do not scale with Q. For Q>,7 GeV, the scaling violations 
are less prononunced, the distributions becoming narrower with increasing Q indicating scaling 
violations of the KNO variable |47j similar to those observed in e + e~ annihilation j|S 



The mean multiplicities, < n c h >, are shown in table [I] as a function of (x, Q 2 ). Over the 
measured range of Q, the mean multiplicity increases by a factor of five. Comparison of the 
mean multiplicity in intervals of fixed Q and differing x provides a check on a possible x 
dependence in the current region of the Breit frame. The mean multiplicities in these intervals 
are the same within systematic errors. 



x range Q 2 range Q mean <n c h> ± stat ± sys 

(GeV 2 ) (GeV) 

0.6 - 1.2 x 10~ 3 10-20 3^8 1.12 ± 0.02 ± 0.15 



1.2-2.4xKT 3 10-20 
20-40 
40-80 

0.24- 1.0 x 10^ 2 20-40 
40-80 
80-160 
160-320 

1.0 - 5.0 x 10~ 2 320-640 
640-1280 



3.8 


1.25 ±0.02 ±0.25 


5.3 


1.73 ±0.04 ±0.23 


7.3 


2.35 ±0.07 ±0.11 


5.3 


1.83 ±0.03 ±0.27 


7.4 


2.47 ±0.05 ±0.23 


10.4 


3.05 ±0.09 ±0.31 


14.5 


3.77 ±0.22 ±0.40 


20.4 


4.49 ±0.30 ±0.54 


29.2 


5.57 ± 0.65 ± 1.21 



Table 1: Mean charged multiplicity in the current fragmentation region. 

In figure 3a the mean charged multiplicity as a function of Q is compared to the predictions 
of three Monte Carlo models which incorporate coherence effects. The MEPS and ARIADNE 
models reproduce the growth in multiplicity as observed in the data. These Monte Carlo 
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generators predict a higher multiplicity for the higher x range at fixed Q values. HERWIG, 
which includes coherent parton showers but has no explicit matrix element, also reproduces 
the increase with Q. Coherence suppresses the production of large multiplicity events which 
becomes a more noticeable effect with increasing Q. 

The effects of coherence on the evolution of the parton shower have been tested by calculating 
the mean multiplicity with the MEPS model (with default parameter settings), as shown in 
figure 3b. The MEPS incoherent parton shower with independent fragmentation exhibits a 
faster growth of the multiplicity than seen in the data. The discontinuities in the lines corre- 
spond to predictions for two different x values at fixed Q. The growth of the incoherent case 
is damped when string fragmentation is used rather than independent fragmentation but not 
enough to agree with the data. Additional low-momentum partons produced in the incoherent 
shower correspond to small kinks of the Lund string. These therefore lead to relatively small 
differences between the incoherent and coherent case when string fragmentation is applied f49 



Overall, however, the data are best described by Monte Carlo models that include coherence 
in their simulations, whilst those with an incoherent shower are less successful in reproducing 
the multiplicity distributions and the growth of <n c h> with Q. 

The definition used here for associating particles with the current fragmentation region, namely 
p z < 0, should be considered an operational definition which is unambiguous only when the 
produced jet and the fragmentation products are massless. If these are massive and/or if QCD 
radiation is included, the current and target fragmentation regions begin to merge. Consider 
for example the case where the struck quark branches into a quark and a gluon producing a (qg) 
system with mass m. Being massive, the ^-momentum of the qg system becomes more positive; 
a special case is m = Q where the qg system is at rest in the Breit frame. If, instead of the 
definition used in this paper, one would take the point of view that the current fragmentation 
should include the particles emerging from the hadronisation of the outgoing quark and gluon 
then the current region defined by p z < would contain only part of the fragmentation products. 
The Bremsstrahlung nature of QCD radiation results in partons which are mostly soft and 
collinear leading to small masses m. Clearly, the effect should be small if Q >> m. Two 
exercises were made in order to estimate the size of the effect, keeping the definition of the 
current region as before, namely p z < 0. 

In the first study, the multiplicities in the current region predicted by the QPM and the ME 
model (i.e. first-order QCD matrix elements without parton showers) were compared. The 
fragmentation was performed in both cases according to the Lund scheme. In the ME calcu- 
lation, QCD branching is considered if y^ > y cut , where = rrijj/W 2 . Here rriy is the 
mass of the parton system produced in the branching and W is the total centre of mass energy 
available for hadron production. The calculation was performed choosing a standard value of 
y C ut = 0.015. Note that this choice of y cut corresponds to a hard branching which is a relatively 
rare process. The average multiplicities predicted by ME at Q = 3.8 GeV were found to be 
~12% lower than those of QPM, the difference becoming smaller as Q increased. 

In the second study, this effect was investigated with data by comparing events with and 
without a large rapidity gap. For large rapidity gap (LRG) events, QCD radiation is suppressed 
compared to events without a large rapidity gap j26fl . Therefore, the depopulation of the current 



region due to QCD radiation should also be suppressed in LRG events. Indeed, in the lowest 
Q interval (Q = 3.8 GeV) LRG events were observed to have an ~ 50% higher multiplicity 
compared with non-LRG; for Q = 7.3 GeV the excess was ~ 15% and for Q > 10 GeV the 
same multiplicity was observed for both event classes, within errors. 
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The preceding considerations suggest that for Q > 10 GeV the different definitions of the 
current region lead to approximately the same results for current fragmentation and that the 
data obtained here from DIS can be directly compared with those on quark fragmentation 
measured in e + e~ annihilation. 

In figure [|, twice the mean multiplicity is compared with inclusive mean multiplicity measure- 
ments from e + e~ experiments |3_8|, [)t], [H], for Q = yfs > 10 GeV. In the region of overlap 
the results from ZEUS and the e + e~ experiments are in agreement, exhibiting a similar rate of 
growth with Q. DIS provides a larger fraction of light quarks compared with the e + e~ case [J3J. 
The different flavour composition of the two physics processes could, in principle, affect the 
comparison of these measurements. However, Monte Carlo generator studies of heavy quark 



production in e + e~ using JETSET and results from e + e _ [|53| on < n c h > from a sample ex- 
cluding 6-quarks show that these effects give rise to differences in the results that are within 
the quoted ZEUS errors. 



8.2 Scaled momentum spectra 

In figure || the \n(l/x p ) distributions for charged particles are shown in the same (x, Q 2 ) intervals 
as the multiplicity distributions. These distributions are approximately Gaussian in shape with 
mean multiplicities given by the integral of the distribution (this estimate agrees within errors 
with the multiplicity results, but differs slightly due to differences in the correction methods). 
The ln(l/x p ) distributions are compared to the ARIADNE Monte Carlo which reasonably 
reproduces the data in each (x, Q 2 ) interval. 

The peak position of the distributions, \ia(l/x p ) max , was evaluated by fitting a Gaussian over 
a range ±1 ln(l/x p ) unit around the statistical mean. This fit was motivated by the MLLA 
prediction for the form of the spectrum (appendix A) which can be approximated by a Gaussian 
distribution around the peak position at sufficiently high energies. For consistency, the same 
fit was performed on TASSO data []TU]] at different centre of mass energies and data from the 
OPAL experiment ULT |. The ln(l/x p ) distributions published in |H| and |ll| include charged 



particles from Kg and A decays: Monte Carlo studies using JETSET indicated that the effect 
of their inclusion on the peak position was less than 0.5%. Results of the fits are shown in 
tables ^| and |||. Our fit to the e + e~ data agrees, within statistical errors, with a similar fit 
performed by the OPAL collaboration [53] Q. 



The ZEUS systematic error includes the contributions from refitting the data with the system- 
atic variations listed in section [7|. These were combined in quadrature with an estimate of the 



uncertainty on the fit obtained by varying the fit range and mean value by up to two bins [5S 
The relative contribution to the systematic error from the fit variation is ~ 40%. 

Figure || shows the distribution of \n(l/x p ) max as a function of Q for the ZEUS data and of yfs 
for the e + e~ data. Over the range shown the peak moves from ~ 1.5 to 2.8, equivalent to the 
position of the maximum of the corresponding momentum spectrum increasing from ~ 400 to 
900 MeV. The ZEUS data complement the e + e~ results. The ZEUS data points are consistent 

lr The OPAL analysis, using a modified Gaussian fit over a limited range around the peak, was not repeated 
due to the limited number of ZEUS data points. A modified Gaussian fit over an extended range leads to 
consistent results for ln(l/x p ) m ax for the ZEUS data |f55| . 
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x range 


Q range 
(GeV 2 ) 


Q mean 
(GeV) 


ln(l/x p ) max ± stat ± sys 


O/i r 

X /dof 


0.6 


- 1.2 x 10" 3 


10-20 


3.8 


1.49 ±0.02 ±0.06 


8.2/6 


1.2 


- 2.4 x 10~ 3 


10-20 


3.8 


1.50 ±0.03 ±0.06 


4.7/6 






20-40 


5.3 


1.70 ±0.03 ±0.06 


14.4/8 






40-80 


7 3 


i 04 + n or + n 1 1 

_L . Jt: _l \J . \J \J _l W . X -L 


4 4/6 


2.4- 


- 10.0 x 10~ 3 


20-40 


5.3 


1.70 ±0.03 ±0.09 


12.5/6 






40-80 


7.4 


1.92 ±0.03 ±0.07 


5.8/6 






80-160 


10.4 


2.18 ±0.06 ±0.10 


3.2/6 






160-320 


14.5 


2.25 ±0.10 ±0.23 


4.0/6 


1.0 


- 5.0 x 10~ 2 


320-640 


20.4 


2.78 ±0.22 ±0.26 


3.8/6 






640-1280 


29.2 


2.85 ±0.21 ±0.24 


5.6/8 



Table 2: Fitted values of \n(l/x p ) max from the ZEUS data. 



Experiment 


yfe (GeV) 


ln(l/x p ) max ± stat 


X 2 /dof 


TASSO 


14 


2.356 ±0.022 


3.9/7 




22 


2.667 ±0.024 


4.1/7 




35 


3.020 ±0.018 


15.3/7 




44 


3.104 ±0.021 


4.6/7 


OPAL 


91 


3.594 ±0.013 


19.1/17 



Table 3: Fitted values of\n(l/xp) ma x from the TASSO and OPAL data. 

with those from TASSO and a clear agreement in the rate of growth of the ZEUS points with 
the e + e~ data at higher Q is observed. 

The increase of ln(l/x p ) max can be approximated phenomenologically by the straight line fit 

\n(l/x p ) max = 6 1n(Q)±c 

also shown in figure |6|. The values obtained from the fit to the ZEUS data are b = 0.650 ± 
0.035(stat) ± 0.069(sys) and c = 0.626 ± 0.059 ± 0.129. The systematic errors are calculated 
by re-fitting the ln(l/a;p) distributions obtained according to the variations listed in section [F] 
and combining the deviations from the central value of the fit parameter, b or c, in quadrature. 
The dominant error is from the correction of the data using the HERWIG Monte Carlo model. 
Removing the lowest two Q points from the fit, where the range over which the Gaussian is fitted 
extends well beyond the momentum range in which the MLLA is valid (see section [5]), leads to 
no statistically significant deviation from the quoted values of b and c. The gradient extracted 
from the OPAL and TASSO data is b = 0.653±0.012 (with c = 0.653±0.047) which is consistent 
with the ZEUS result. This value is consistent with that published by OPAL, b = 0.637 ±0.016, 
where the peak position was extracted using an alternative method |TT|]. A consistent value of 



the gradient is therefore determined in DIS and e e annihilation experiments. 
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Also shown is the statistical fit to the data when b = 1 (c = 0.054 ± 0.012) which would 
be the case if the QCD cascade was of an incoherent nature, dominated by cylindrical phase 
space. (A discussion of phase space effects is given in [jnj.) In such a case, the logarithmic 
particle momentum spectrum would be peaked at a constant value of momentum, independent 
of Q. The observed gradient is clearly inconsistent with b = 1 and therefore inconsistent with 
cylindrical phase space. 

In figure [F|a the measurements are compared with the results from the ARIADNE, HERWIG 
and MEPS Monte Carlo simulations. These models describe the \n(l/x p ) distributions and the 
evolution of their maximum as well as the multiplicity distributions shown in section ETT1 . As in 
the case of < n c h >, the difference between the coherent and incoherent growth of hx(l/x p ) max is 
found to be reduced when string fragmentation is used rather than independent fragmenta- 
tion (figure |7|b). Again, the ln(l/x p ) distributions and the growth of \n.(l/x p ) max with Q are 
well- described by Monte Carlo models that include coherence in their simulations, whilst those 
with an incoherent shower are less successful in describing the data. 

In contrast to the < n c h > measurements, the \n(l/x p ) max values change only slightly when only 
the rapidity gap events are used, with a maximum negative shift of ~ 5% in the lowest (x, Q 2 ) 
interval. Effects due to primary heavy-quark production were also found to be negligible. 

The evolution of ln(l/ 'x p ) max seen in DIS agrees with the e + e~ data. The growth is at a slower 
rate than that expected from cylindrical phase space, indicative of coherence. Monte Carlo 
models including coherence effects describe the ln(l / x p ) max values and the ln(l/x p ) distributions 
better than those with an incoherent behaviour. The rate of growth with increasing Q of < n c h > 
and ln(l/x p ) max in the current region of the Breit frame agrees with e + e~ data in the region 
of overlap (i.e. above 10 GeV) indicating that the fragmentation of quarks in DIS is similar to 
that of quarks produced in e + e~ annihilation. 

9 Comparisons with MLLA 

In order to test the predictions of the MLLA (Modified Leading Log Approximation) and LPHD 
(Local Parton Hadron Duality) at different energies, fits to the individual \n(l/x p ) distributions 
were performed in the chosen (x, Q 2 ) intervals. The evolution of these distributions with Q is 
also predicted within the framework of the MLLA. 

In the high energy limit, MLLA calculations predict the momentum distribution of soft gluons, 
D hm , radiated by a quark of energy E = Q/2 according to equation (f|) given in appendix A. 
D hm is known as the limiting spectrum. Invoking LPHD, the hx(l/x p ) distribution of the 
hadrons is given by 



where Y = ln(Q/2A). The A in the expression is an effective scale parameter and not the 
standard QCD scale, e.g. Aj^r. This MLLA prediction is valid in the range 0<,ln(l/x p ) <Y, a 
region where pQCD can be applied. The normalisation factor K ch is energy independent within 
the framework of LPHD. In addition, the charged multiplicity can be expressed as 



1 da 



K ch D lim (\n(l/xp),Y) 



(1) 



a d(ln(l/xp)) 



<n ch >= K ch Af(Y) 



(2) 
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where M{Y) is the integral of D hm given by equation (0) in appendix A and represents the soft 
gluon multiplicity. 

The evolution of \n(l/xp) max with Q is a function of A and is given by 

Hl/x p ) max = ^Y + c 2 W -cj+ (OiY- 1 / 2 )) (3) 

where c 2 is a constant calculated from the number of flavours and colours. Taking the active 
number of flavours to be 3, the value of the constant c 2 is 0.52. The term 0(Y~ 1 / 2 ) describes 
higher-order effects. 

The functional form of equation ([!]) was fitted to the ln(l/x p ) distributions. This equation has 
two free parameters, n ch and A. The fit was performed for each (x, Q 2 ) interval in the range 
< ln(l/x p ) < \n(Q/2p ) where p = 400 MeV. This restriction limits the comparisons of 
MLLA at the lower Q values. Due to the restricted range of the fit very little of the peak region 
of the lowest Q value is included: these intervals are therefore excluded from the following fits. 
The results of these fits are shown in table [|. 





x range 


Q 2 range 
(GeV 2 ) 


Q mean 
(GeV) 


A± stat 
(MeV) 


n ch ± stat 


1.2 


-2.4 x 10~ 3 


20-40 


5.3 


330 ± 29 


1.25 ±0.05 






40-80 


7.3 


323 ± 20 


1.28 ±0.07 


2.4- 


- 10.0 x 10~ 3 


20-40 


5.3 


313 ±28 


1.21 ±0.04 






40-80 


7.4 


296 ± 19 


1.30 ±0.05 






80-160 


10.4 


296 ± 28 


1.25 ±0.06 






160-320 


14.5 


311 ±48 


1.31 ±0.10 


1.0 


- 5.0 x 10~ 2 


320-640 


20.4 


210 ±54 


1.14 ±0.10 






640-1280 


29.2 


266 ±81 


1.18 ±0.16 



Table 4: Fitted values of A and k' 



Approximately the same value of A is found at all values of Q. The values of A are combined 
to give a weighted mean value of 

A = 306 ± 10 (stat) ± 30 (sys) MeV. 

The fits in each interval are also consistent with a single value of K ch with a mean value of n ch 
given by 

K ch = 1.25 ± 0.02 (stat) ± 0.09 (sys). 

The experimental systematic errors are calculated by re-fitting the ln(l/x p ) distributions ob- 
tained according to the variations listed in section [I] and combining the deviations from the 
central values of the fit parameters, A and K ch , in quadrature. The dominant error is from 
the correction of the data using the HERWIG Monte Carlo model. The absolute value of the 
correlation coefficient between A and n ch from the fits was typically ~ 0.35. Equation ([!]) was 
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re-fitted with p varying from 200 to 500 MeV, yielding consistent results within the quoted 
systematic errors. Similarly, this equation was re-fitted with A fixed at its mean value and K ch 
as a free parameter. The values of n ch remain the same within the statistical accuracy of the 
fit with little degradation in the \ 2 . 

The results of the fits to the MLLA limiting spectrum using the mean values of A and n ch 
above are shown in figure |[ The data are reasonably well described by the limiting spectrum 
with single values of A and K ch over a large range in Q. The overall x 2 /dof for all the ZEUS 
data is 152/106. At higher values of ln(l/x p ), beyond the range of the fit, the MLLA spectrum 
underestimates the multiplicity. 

The values measured by OPAL, fitting their data over a limited momentum range, were 
A = 253 ± 30 MeV and K ch = 1.44 ± 0.01. In order to compare with the ZEUS measurement, 
their value of K ch has been divided by a colour factor of 4/9 and by a factor of 2 to account for 



the quark and anti-quark pair, factors which OPAL absorbed into their value quoted in [II 
The inclusion of decay-particles from K$ and A in the OPAL data contributes 10% to the mul- 



tiplicity |jT| and can account for the higher value of k . The value of A is also sensitive to this 
but to a lesser degree. Monte Carlo studies using JETSET indicated that the effect of their 
inclusion on the value of A was less than 10 MeV. 

In order to check the consistency of K ch , equation (|2]) was fitted to < n c h > as a function of 
Q. Taking the value of A extracted from the individual ln(l/x p ) distributions, A = 306 MeV, 
the fitted value of K ch = 1.78 ± 0.02 (stat) was determined. This value of n ch is significantly 
higher than that extracted from the ln(l/x p ) distributions. Equation (|2|) is derived from the 
integration of equation (D), this expression underestimates the data at large ln(l/x p ) (see 
figure |8|). At the relatively low Q values of the ZEUS data, Q^7 GeV, this accounts for a 
significant fraction of the total multiplicity where the MLLA prediction, taken outside its range 
of validity, underestimates the data. The underestimate of the predicted parton multiplicity by 
the MLLA is therefore compensated by an increase in the value of K ch . 

In order to check the consistency of A, the evolution of \n(l/x p ) max was investigated. The 
value of A hardly influences the slope of \n(l/x p ) max in equation but gives the absolute 
value of ln(l/x p ) max at a given Q value. This is highly correlated with the higher-order term 
(9(V -1 / 2 ) in equation (W) and therefore this term is neglected in the following fits. The Gaussian 



fits which determine the \n(l/x p ) max values described in section |8.2| extend to momentum 
values lower than the measured A. In order to compare with MLLA calculations, the peak 
position of the ln(l/x p ) distribution was extracted again, using a Gaussian symmetric around 
the estimated peak position. The range of the fit was from (2L — Y) to Y where L is 
the estimated peak position, calculated from equation (|3|) neglecting the 0(Y~ 1 ^ 2 ) term and 
assuming Y = \n(Q/2p ), where po = 306 MeV. A value of A = 284 ± 8(stat) MeV is found, 
consistent with the value extracted from the fit of the individual plots. This compares to a 
value of A = 212 ± 20 MeV from the evolution of the TASSO and OPAL data as extracted 
by OPAL. (OPAL fit equation (|3p with c<i as a free parameter; a high correlation was found 
between C2 and A when fitting the ZEUS data.) 

Although the growth of the multiplicity cannot be consistently described by the MLLA and 
LPHD, a single value of A and K ch can describe the ln(l/;Cp) distributions over a range of 
Q values from ~ 5 to 30 GeV. The values obtained are consistent with those extracted at 
y/s = 91 GeV by the OPAL experiment. The evolution of ln(l/x p ) max with Q also gives a 
consistent value of A, again compatible with the OPAL and TASSO results. 
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10 Conclusions 



Charged particle distributions have been measured in the current region of the Breit frame 
in DIS over a wide range of Q values. In the framework of the Monte Carlo models, the best 
description of the data is achieved by those incorporating coherence effects where the production 
of soft particles and the growth of the mean multiplicity is suppressed relative to the incoherent 
case. Modified Leading Log Approximation (MLLA) calculations have been shown to describe 
the \n(l/xp) distributions in the range 5 <, Q < 30 GeV using single values of n ch , which relates 
the observed hadron distributions to the calculated parton distributions, and of A, an effective 
scale parameter of the QCD calculations. 

A comparison of these charged particle distributions with those from e + e~ annihilation exper- 
iments has been performed. For Q > 10 GeV a growth in <n c h> is observed with Q which is 
similar to the growth in e + e~ experiments as a function of y/s. The evolution of ln(l/x p ) max has 
been measured and is consistent with that observed in e + e~ experiments. The fragmentation 
properties of the struck quark from the proton in DIS that have been studied here are similar 
to those from quarks created in e + e~ annihilation. The observed charged particle spectra are 
therefore consistent with the universality of quark fragmentation. 
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Appendix A 
MLLA formulae 

QCD predicts the form of the soft gluon momentum distribution at a fixed energy and describes 
the energy evolution of the spectrum, including the multiplicity and the peak maximum. 

In the high energy limit, MLLA (pQCD) calculations predict || the momentum distribution 
of soft gluons from a quark with energy E according to the following formula, in which xp is 
the momentum scaled with E, 



where C/, b and B are constants based on the number of flavours and colours. For the number 
of flavours, Nf = 3, and the number of colours, N c = 3, the constants Cf, b and B are equal to 
4/3, 9, and 1.247 respectively. Y = \a(E/A), where A is an effective scale parameter. For the 
case of DIS, in the current fragmentation region, E = Q/2 and for e + e~, E = y/s/2. F denotes 
the Gamma function and Ip is a modified Bessel function of order B. The variable a is given 
by a = cto + ir, where ao is determined by 



and ( = 1— ln(l/xp)/Y. The integral, which is a representation of the confluent hypergeometric 
function, is performed over r. D hm is known as the limiting spectrum which is expressed in a 
form that is suitable for numerical integration. 

From the limiting spectrum the soft gluon multiplicity can be calculated as 




1 B /2 




tanh cto = 2£ — 1 




(5) 



where z 



b 
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Figure 1: Charged multiplicity distributions in the current region of the Breit frame as a function 
of (x, Q 2 ). Statistical errors are indicated by the inner error bar bounded by the horizontal bars. 
The outer error bars show the statistical and systematic errors added in quadrature. The ZEUS 
data are compared to the predictions of the ARIADNE Monte Carlo. 
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Figure 2: TTie ifiVO distribution in the current region of the Breit frame as measured by ZEUS. 
The values of Q correspond to the x-intervals shown in table 1. The highest (x,Q 2 ) interval is 
omitted for the sake of clarity. Statistical errors only are indicated. 
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Figure 3: <n c / l > as a function of Q for the x-intervals shown in table 1. The ZEUS data are 
compared to various Monte Carlo models. The data in figures (a) and (b) are the same. The 
discontinuities in the Monte Carlo curves correspond to predictions for two different x values 
at fixed Q. The line connects the lower-x to the higher-x point with increasing Q. 
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Figure 4: Mean charged multiplicity as a function ofQ. Twice the measured ZEUS multiplicity, 
2- < n ch >, is compared to < n ch > results from HRS, OPAL, PLUTO and TASSO. All data 
were corrected for K s and A decays. 
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ZEUS 1993 




Figure 5: \n(l/x p ) distributions in the current region of the Breit frame as a function of(x,Q 2 ). 
Statistical errors are indicated by the inner error bar bounded by the horizontal bars. The outer 
error bars show the statistical and systematic errors added in quadrature. The ZEUS data are 
compared to the predictions of the ARIADNE Monte Carlo. 
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Figure 6: ln(l/x p ) max as a function of Q for the x-intervals shown in table 1. The 
ZEUS data are compared to results from OPAL and TASSO. A straight line fit of the form 
\n(l/x p ) max = b ln(Q) + c to the ZEUS \n(l/x p ) max values is indicated as well as the line 
corresponding to b = 1, discussed in the text. 
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Figure 7: ln(l/x p ) maa; as a function ofQ. The ZEUS data are compared to various Monte Carlo 
models. The data in figures (a) and (b) are the same. The discontinuities in the Monte Carlo 
curves correspond to predictions for two different x values at fixed Q. 
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Figure 8: ln(l/x p ) distributions in the current region of the Breit frame as a function ofQ 2 . The 
ZEUS data are compared to the MLLA limiting spectrum in the range < ln(l/x p ) < ln(Q/2A) 
with A = 306MeV and n ch = 1.25. For clarity, only the lower x range is shown in the two lowest 
Q 2 intervals. 
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